fixedSubs2 <- function(chick.transcript,grouse.transcript,dir="/Paterson/Datafiles/grouse/test_files",transition.ratio=FALSE,sliding.window=10,max.diffs=6,trim.seq=FALSE){
	#calculates fixed differences between a pair of sequences
	#(called chicken and grouse)
	require(GeneR)
	muscle.char <- character(4)
	muscle.char[1] <- paste(">ChickenProtein")
	muscle.char[2] <- strTranslate(chick.transcript)
	muscle.char[3] <- paste(">GrouseProtein")
	muscle.char[4] <- strTranslate(grouse.transcript)
	
	muscle.file <- paste(dir,"/proteins_tmpfile.fa",sep="")
	muscle.file.out <- paste(dir,"/proteins_tmpfile.afa",sep="")
	writeLines(muscle.char,muscle.file)
	
	ta.char <- character(4)
	ta.char[1] <- paste(">ChickenTranscript")
	ta.char[2] <- chick.transcript
	ta.char[3] <- paste(">GrouseTranscript")
	ta.char[4] <- grouse.transcript
	
	tranalign.file <- paste(dir,"/transcripts_tmpfile.fa",sep="")
	tranalign.file.out <- paste(dir,"/transcripts_tmpfile.afa",sep="")
	writeLines(ta.char,tranalign.file)

	#need to set R up use same PATH as X11
	system(paste("/usr/local/bioinf/muscle/muscle -in",muscle.file,"-out",muscle.file.out))
	system(paste("/opt/local/bin/tranalign",tranalign.file,muscle.file.out,tranalign.file.out))
	require(Biostrings)
	prot.afa <- read.BStringSet(muscle.file.out)
	dna.afa <- read.BStringSet(tranalign.file.out)
	
	
	tst <- strsplit(as.character(prot.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	tst2 <- list(chick.prot = toupper(tst[[1]][-rm.gaps]), grouse.prot = toupper(tst[[2]][-rm.gaps]))
	
	#insert homology criteria here
	#
	prot.homology <- which(tst2$chick.prot==tst2$grouse.prot)
	prot.nonhomo <- which(tst2$chick.prot!=tst2$grouse.prot)
	prot.exclude <- rep(FALSE,length(tst2$grouse.prot))
	
	if(prot.homology[1]>1){
		prot.exclude[1:(prot.homology[1]-1)] <- TRUE
		}
	if(prot.homology[length(prot.homology)]<length(tst2$grouse.prot)){
		prot.exclude[length(tst2$grouse.prot):(prot.homology[length(prot.homology)]+1)] <- TRUE
		}
	
	#sliding window for differences within alignment
	#sliding.window <- 10 #change this if needed
	#max.diffs <- 6 #change this if needed
	diff.list <- vector('list',length=length(tst2$grouse.prot)-sliding.window+1)
	for(ipos in 1:length(diff.list)){
		diff.list[[ipos]] <- prot.nonhomo[prot.nonhomo >=ipos & prot.nonhomo<ipos+sliding.window]
		}
	diff.count <- sapply(diff.list,length)
	diff.unlisted <- unique(unlist(diff.list[diff.count>max.diffs]))
	#diff.unlisted should contain vector of sites to exclude
	#note, this will only exclude non-homologous sites, not homologous sites
	#embedded within a group of non-homologous sites
	prot.exclude[diff.unlisted] <- TRUE
	
	if(trim.seq) {
		tst2$chick.prot <- tst2$chick.prot[!prot.exclude]
		tst2$grouse.prot <- tst2$grouse.prot[!prot.exclude]
		dna.char <- strsplit(as.character(dna.afa),split="")
		rm.gaps <- unique(c(grep("-",dna.char[[1]]),grep("-",dna.char[[2]])))
		if(length(rm.gaps) > 0){
			dna.char[[1]] <- dna.char[[1]][-rm.gaps]
			dna.char[[2]] <- dna.char[[2]][-rm.gaps]
			}
		
		wpe <- which(prot.exclude)
		if(length(wpe) > 0){
			dna.exclude <- c((3*wpe)-2, (3*wpe)-1, 3*wpe)
			dna.char[[1]] <- dna.char[[1]][-dna.exclude]
			dna.char[[2]] <- dna.char[[2]][-dna.exclude]
			}
		dna.names <- names(dna.afa)
		dna.afa <- DNAStringSet(c(paste(dna.char[[1]],collapse=""),paste(dna.char[[2]],collapse="")))
		names(dna.afa) <- dna.names
		}
	
	
	#yn00 estimates
	paml.estimates <- try(pamlCalc(dna.afa,dir))
	#note, if paml fails may still take results from existing files
	
	prot.diff <- sum(tst2$chick.prot!=tst2$grouse.prot)
	tst <- strsplit(as.character(dna.afa),"")
	rm.gaps <- unique(c(grep("-",tst[[1]]),grep("-",tst[[2]])))
	if(length(rm.gaps)>0){
		tst2 <- list(chick.dna = toupper(tst[[1]][-rm.gaps]), grouse.dna = toupper(tst[[2]][-rm.gaps]))
		}else{
		tst2 <- tst #this is getting to be a mess of patches
		names(tst2) <- c("chick.dna","grouse.dna")
		}
	aligned.length <- length(tst2$grouse.dna)
	
	if(any(!tst2$chick.dna %in% c("A","G","C","T"))) warning("unexpected base in chicken sequence\n")
	if(any(!tst2$grouse.dna %in% c("A","G","C","T"))) warning("unexpected base in grouse sequence\n")
	dna.diff <- sum(tst2$chick.dna!=tst2$grouse.dna)
	if(transition.ratio==TRUE){
		#count up number of transitions and tranversions
		tst2$chick.type <- sapply(tst2$chick.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in chicken sequence')
			return("N") #only if function gets this far
			})
		tst2$grouse.type <- sapply(tst2$grouse.dna,function(X){
			if(X %in% c("A","G")) return("R") #purine
			if(X %in% c("C","T")) return("Y") #pyrimidine
			warning('unknown base in grouse sequence')
			return("N") #only if function gets this far
			})
		no.transversions <- sum(tst2$chick.type!=tst2$grouse.type)
		no.transitions <- dna.diff - no.transversions
		tmp <- c(aligned.length,prot.diff,dna.diff-prot.diff,no.transitions,no.transversions)
		names(tmp) <- c("aligned.length","Dn","Ds","transitions","transversions")
		}else{
		tmp <- c(prot.diff,dna.diff-prot.diff)
		names(tmp) <- c("Dn","Ds")
		}
	if(class(paml.estimates)=="try-error"){
		tmp <- c(tmp,rep(-1,9))
		}else{
		tmp <- c(tmp,unlist(paml.estimates))
		}
	tmp
	}

grouse.chick.MK2 <- grouse.chick.map
grouse.chick.MK2[,c("aligned.length","Dn", "Ds", "transitions", "transversions", "dN", "dN.se", "dS", "dS.se", "kappa", "omega", "Ls", "Ln", "t")] <- -1
#Sys.setenv(PATH=paste(Sys.getenv("PATH"),"/usr/local/bioinf/paml/bin",sep=":"))


for(grouse.gene in 1:nrow(grouse.chick.map)){#change to nrow(grouse.chick.map)
	chick.transcript <- as.character(chick.seqs[[grep(grouse.chick.map$ensembl_transcript_id[grouse.gene],names(chick.seqs))]])
	grouse.transcript <- getCDS(seq.name=grouse.chick.map$grouse.id[grouse.gene],single.cds=TRUE,dir.embl="/Paterson/Datafiles/grouse/embl_exon3")
	if(!"transcript" %in% names(grouse.transcript)) next
	
	#remove *s from protein/transcript
	star.pos <- grep('\\*',strsplit(grouse.transcript$protein,'')[[1]])
	if(length(star.pos)>0){
		tst.trans <- strsplit(grouse.transcript$transcript[[1]],'')[[1]]
		cut.trans <- tst.trans[-c( (star.pos-1)*3+1, (star.pos-1)*3+2, (star.pos-1)*3+3 )]
		grouse.transcript$transcript <- paste(cut.trans,collapse="")
		}
	
	test.MK <- try(fixedSubs2(chick.transcript,grouse.transcript$transcript,transition.ratio=TRUE,trim.seq=TRUE))
	if(class(test.MK)=="try-error") next
	grouse.chick.MK2[grouse.gene,c("aligned.length","Dn", "Ds", "transitions", "transversions", "dN", "dN.se", "dS", "dS.se", "kappa", "omega", "Ls", "Ln", "t")] <- test.MK
	}

#add polymorphisms
grouse.chick.MK2$Pn <- -1
grouse.chick.MK2$Ps <- -1
grouse.chick.MK2$fisher.pv <- -1
grouse.chick.MK2$fisher.or <- -1
grouse.chick.MK2$fisher.lower <- -1
grouse.chick.MK2$fisher.upper <- -1

for(mki in 1:nrow(grouse.chick.MK2)){
	if(!grouse.chick.MK2$grouse.id[mki] %in% names(grouseSNPdfs)) next
	tmp.snp <- QCsnp(grouseSNPdfs[[grouse.chick.MK2$grouse.id[mki]]])
	if(nrow(tmp.snp)==0){
		grouse.chick.MK2[mki,c("Pn","Ps")] <- c(0,0)		}else{
		tmp.cnt <- colSums(tmp.snp[,c("cds","dn")])
		grouse.chick.MK2[mki,c("Pn","Ps")] <- c(tmp.cnt['dn'],tmp.cnt['cds']-tmp.cnt['dn'])
		tmp.fisher <- try(fisher.test(matrix(as.numeric(grouse.chick.MK2[mki,c('Dn','Ds','Pn','Ps')]),nrow=2,byrow=TRUE)))
		if(class(tmp.fisher)=="try-error"){
			cat(mki,"  ",grouse.chick.MK2$grouse.id[mki],'  failed')
			next
			}
		grouse.chick.MK2$fisher.pv[mki] <- tmp.fisher$p.value
		grouse.chick.MK2$fisher.or[mki] <- tmp.fisher$estimate
		grouse.chick.MK2$fisher.lower[mki] <- tmp.fisher$conf.int[1]
		grouse.chick.MK2$fisher.upper[mki] <- tmp.fisher$conf.int[2]
		}
	
	}

#pass data to Welch test
#columns are:
#1   2     3     4       5   6      7   8     9
#Dn	LN(D) PN	LN(P)	DS	LS(D)	PS	LS(P) alleles

#note that LN/S(P) will be a slight under-estimate, some gapped sites 
#in alignment may be excluded. But sites present in grouse but not
#chicken tend to be rare. LN/S(P) taken from untrimmed alignments.

grouse.chick.MK[grep('Error',grouse.chick.MK$Ls),"Ls"] <- "-1"
grouse.chick.MK[grep('Error',grouse.chick.MK$Ln),"Ln"] <- "-1"
grouse.chick.MK[grep('Error',grouse.chick.MK$Ln),"omega"] <- "-1"

grouse.chick.MK$Ls <- as.numeric(grouse.chick.MK$Ls)
grouse.chick.MK$Ln <- as.numeric(grouse.chick.MK$Ln)

Welch.trimmed <- cbind(grouse.chick.MK2[,c("grouse.id","ensembl_transcript_id","Dn","Ln","Pn")],
 	grouse.chick.MK[,"Ln"],grouse.chick.MK2[,c("Ds","Ls","Ps")],grouse.chick.MK[,"Ls"])
names(Welch.trimmed) <- c("grouse.id","ensembl_transcript_id","Dn","LnD","Pn","LnP","Ds","LsD","Ps","LsP")


tmp <- data.frame(cDNA=grouse.df3$cDNA , depth=rowSums(grouse.df3[,grep('depth',names(grouse.df3))]),stringsAsFactors=FALSE)

Welch.trimmed$alleles <- 100 #should take this from tsv coverage?

Welch.trimmed2 <- Welch.trimmed[Welch.trimmed$Dn!= -1,]

write.table(Welch.trimmed2[,3:ncol(Welch.trimmed2)],file='/Paterson/Datafiles/grouse/test_files/Welch_test1.csv',quote=FALSE,row.names=FALSE,col.names=FALSE,sep=",")